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We study a single impurity driven by a constant force through a one-dimensional Bose gas using a Lieb- 
Liniger based approach. Our calculaton is exact in the interaction amongst the particles in the Bose gas, and is 
perturbative in the interaction between the gas and the impurity. In contrast to previous studies of this problem, 
we are able to handle arbitrary interaction strength for the Bose gas. We find very good agreement with recent 
experiments [Phys. Rev. Lett. 103 , 150601 (2009)]. 


I. INTRODUCTION 

The study of quantum mechanical systems out of equilib¬ 
rium is one of the great frontiers of modern physics. The 
questions in this field are not only of fundamental interest, 
but are also of interest to future quantum technologies, as well 
as to classical technologies on the nano-scale. Cold atomic 
systems have provided an ideal setting for hand-in-hand the¬ 
oretical and experimental investigations of this frontier, par¬ 
ticularly in low dimensions. Nonetheless, our understanding 
of the issues involved are sufficiently primitive that it remains 
useful to consider some of the simplest toy model experiments 
in order to gain intuition regarding more general questions. 

In this paper we focus on the study of driven impurities 
moving through a one-dimensional (ID) Bose gas. This sub¬ 
ject has received much attention of late, thanks to both ex¬ 
perimental and theoretical progress |@-[2^ Our work 
was inspired by the experimental results in Ref. [J, where the 
impurity is driven through the gas by a constant force (grav¬ 
ity). This type of experiment has not been extensively investi¬ 
gated at the theoretical level, with the exception of the recent 
work in Re f. [H (see also Ref. [20h which uses a Tonks-Gas 
description 0^ (appropriate only in the limit of large inter¬ 
action strength). We also refer the reader to Ref. HI where a 
similar system was studied in presence of a ID optical lattice. 

In this work, we use linear response theory and Fermi’s 
golden rule with exact transition rates to model the scatter¬ 
ing between the driven impurity and the underlying gas (see 
e.g., Refs. d[Hill. We then model the motion of the impu¬ 
rity as a classical driven stochastic process lIH . Our approach 
is strictly valid in the limit where the interaction between the 
impurity and the underlying gas is sufficiently weak (which is 
not necessarily true in the experiments of Ref. Qj)- However, 
in contrast to prior works attempting to analyze this problem, 
our approach is valid for any interaction strength between par¬ 
ticles in the ID gas. The main point of this work is to provide a 
method to analyze the effects of interaction within the ID gas 
on a driven impurity - which has previously not been possible 
for intermediate interaction strength within the gas. 

We find quantitative agreement with the results in Ref. Q] 
both comparing the center of mass motion as well as the opti¬ 
cal density profile of a packet of impurities after the\^ave the 
ID gas. In contrast to earlier theoretical work lU^I^ which 


found agreement in the strong interaction limit with Tonks- 
Girardeau modelling, our results quantitatively describe the 
experimental data for small values of the dimensionless inter¬ 
action strength 7^1 — 3, close to the Bogoliubov limit liHl . 
in keeping with the experimental parameter range (see Ta- 
blelD. 


II. METHOD 


Our general method for handling driven impurtity motion is 
to treat the scattering of the impurity atom perturbatively. Our 
second key assumption is that the driven impurity atom is a 
negligible perturbation of the underlying ID system, which is 
assumed to relax to its ground state between any two scatter¬ 
ing events. This assumption is strictly valid when the impurity 
scatters only once in the entire time span of the experiment 
(which is approximately the case if the coupling of the im¬ 
purity to the ID gas is weak). Nonetheless, there are several 
other regimes for which this approximation is expected to be 
quite good. For example, the behaviour of the ID gas may 
not differ much if it is in its exact ground state versus being 
slightly excited. Another regime of interest is when the im¬ 
purtity moves faster than the effective speed of propagation in 
the ID gas. In this case if the impurity scatters a second time, 
it will have out-run the perturbation it caused in the first scat¬ 
tering event and will effectively see the ID gas as if it were in 
its ground state. 

We consider a delta function interaction potential of inter¬ 
action strength gim, V = ~ ^0’ between the 

driven impurity atom at position x and the ID gas atoms at 
positions {xi}. In the usual way, Fermi’s golden rule gives a 
transition rate between an initial |i) and final |/) state of the 
system, of energies and Ej respectively, as 
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where the superscript ° indicates that these states and ener¬ 
gies are to be evaluated in the absence of the coupling V be¬ 
tween the impurity and the gas. Hence we have |j°), |/°) = 
|fc) (g) |n), where iFiDgas|n) = e„|n) and iFimpurity|fc) = 
/2mim)\k) with rriim the impurity mass, iTimpurity the 
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Hamiltonian of the impurity, -HiDgas and e„ the Hamiltonian 
and eigenenergies of the ID gas. 

As described above, we assume that only the ground state 
n = 0 appears in |i°). Summing over all final states of the ID 
gas, we obtain a transition rate for the impurity 
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where pk'-k = is the Fourier transform of the 

density operator of the ID gas and S{k, to) is the dimension¬ 
less dynamic structure factor (DSF) of the ID gas (note the 
factor of ep/Np in our definition of S, Np being the number 
of particles and ep = h^kp/2m their Fermi energy; here m 
is the mass of the particles in the ID gas and kp = irn is the 
Fermi wavevector with n = Np/L the density). 

We assume our ID gas is made of spinless bosons and has 
short ranged interactions of the form gip J2i<j ~ ^j)- 
For convenience, we introduce the standard dimensionless in¬ 
teraction strength 7 = gipm/{h^n). Analytical solutions 
for S{k,uj) are available in the weakly or strongly interact¬ 
ing limit. For intermediate values of 7 one may use the exact 
Lieb-Liniger (LL) solution for the DSF which can be obtained 
numerically for any values of k and w. A description of this 
numerical procedure can be found in Ref.HH 

Once we can calculate the transition rate, we need to ac¬ 
count for the driven motion of the impurity. To simulate 
both the driving force and the scattering, we discretize time 
and momentum and write a scattering transition probability 
Vk^k' = Wk^k'StAk and we define the probability of not 
scattering to be Vk^k = 1 - Y^k’^k 'Pk^k'- For every time 
interval 5t, we evolve the position and velocity of the particles 
deterministically. In the present case (inspired by the exper¬ 
iments of Ref. [2) we are concerned with the impurity being 
accelerated (driven) by gravity g (assume acceleration in the 
-fa; direction) so we have 


parameters 

(1) 

(2) 

(3) 

Np in central ID system 

36 

32 

30 

L of central ID system (pm) 

22.46 

24.1 

24.82 

7ave average over entire condensate 

3 

5 

7 

7ct average over central ID system 

1.7 

2.9 

4.0 

7 c at the center of condensate 

1.1 

1.9 

2.6 


TABLE I. Range of parameters of the ID systems obtained from 
three different Bose condensates (columns 1-3) fISl . 


system is in an effective zero-temperature ID regime. The pa¬ 
rameters of the ID traps vary with position - both between 
different ID systems in the ensemble and within each indi¬ 
vidual ID system. It should be noted that the value of 7 is 
expected to vary significantly across the system. Private com¬ 
munication with the authors of Ref. [I] lead to the estimates 
reported in Table |I] 

For simplicity we crudely neglect the nonuniformity of the 
system, considering only the case of a homogeneous ID sys¬ 
tem with fixed density IXh . 

A radio-frequency (RF) pulse is used to change the hy- 
perfine ground state of some (up to 3) atoms near the cen¬ 
ter of the ID system so as to decouple them from the trap. 
The pulse is Fourier-limited in width [full width at half max¬ 
imum (FWHM) ^ 2.3 pm] and has a velocity distribution 
of width ~ 2 • 10 “^ m/s (close to the uncertainty limit). 
Therefore, we consider a wave packet '0(a;) cx sin (ax) jx, 
where a = 27r/Ax ~ 2.73 pm“^ (although we find that 
a ~ 1.8 — 2.0 pm“^ produces a better fit to the width of 
the unscattered peak in the experimental density profile of the 
falling atoms at long times, illustrated in Fig.|2l). 

We model these initial conditions using a Gaussian- 
smoothed ll28l] Wigner function for the position and mo¬ 
mentum distribution of the falling atoms at time t = 0 : 
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x{tdt) = x{t)v{t)6tgSt^/2 
v{t -f 6t) = v{t) -b gSt. 

After each time interval 6t we then allow for a stochastic scat¬ 
tering attempt vi = v{t 6t) —>■ V 2 with probability per unit 
wave vector Vvi^v 2 = Wk^k' dt, where vi = hk/rriim and 
V 2 = hk'/rriim- This allows for efficient simulation of the im¬ 
purity motion. In the large and small 7 regime, we have used 
analytic forms of the DSF to test our numerical algorithm. 


III. EXPERIMENTAL PARAMETERS 

As demonstration of our method we apply it to the exper¬ 
imental situation from Ref. Ql A Bose condensate of ®^Rb 
atoms is confined into an ensemble of harmonic traps with 
long axis aligned with the Earth gravitational field. The trans¬ 
verse radius of each trap and temperature are such that each 


where a and [3 are positive real constants that satisfy the con¬ 
dition ajd < 1, i.e., the smoothing area is > h. We choose the 
least possible smoothing that yields a positive semi-definite 
probability, namely a/3 = 1. The value of a is then set to 
equal ha^. After a few lines of algebra one obtains: 
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where Erf is the Gaussian error function extended to the com¬ 
plex plane. This distribution was sampled using the rejection 
sampling technique. 

The decoupled impurity atoms are then allowed to acceler¬ 
ate under the constant driving force of the gravitational field. 
In our simulations we assume only a single impurity atom is 
decoupled from the trap (i.e, we neglect interactions between 
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multiple falling impurity atoms; we also disregard possible ef¬ 
fective interactions between the decoupled atoms that may be 
mediated by the condensate). 

In this experiment, the falling (impurity) atoms are identical 
to the atoms in the trap (gas) up to their spin state. Hence, 
'm-im = fn and all interactions (impurity-gas and gas-gas) are 
described by the same delta function potential, gim = gin- 


IV. CENTER OE MASS 

We start by considering the position of the overall center 
of mass of the packet of falling atoms as a function of time, 
which was measured experimentally and reported in Fig. 3 of 
Ref.[l]. The parameters used in this experiment are those listed 
as case (3) in Table |I] 

Fig Ushows a comparison between the experimental curves 
from Ref. [l] and the results from our stochatic simulations, 
using n = Np/L ~ 1.2 and different values of 7 (top 

panel), as well as using 7 = 2.6 and different values of n 
(bottom panel). In the simulations we consider an infinite ID 
gas of uniform density. The experimental time f = 0 in Ref.Q] 
was chosen to correspond to the middle of the RF pulse that 
creates the packet of falling atoms. Accordingly, we chose 
t = 0 in our simulations as the time when the Fourier limited 
packet starts moving through the ID gas. 

The numerical results appear to be very sensitive to the val¬ 
ues of 7 and n. This allows us to determine that the combi¬ 
nation n ~ 1.2 fim~^ and 7 = 2.6 provides the best fit to 
the experimental data. (We note that particles begin to leave 
the ID trapped gas after about 2 ms, which corresponds to 
the longest time reported in Fig. [T] Such an effect may be 
responsible for the discrepancy that we observe between our 
numerics and the very last data point in the experiment.) 

Our results are in contrast with earlier theoretical mod¬ 
elling iQl, [Isl [ 2 ^ which achieved a similarly good fit to the 
experimental results by using the strongly interacting Tonks- 
Girardeau (TG) approximation corresponding to 7 = 00 
within the ID gas and then treating the interaction between 
the impurities and the ID gas at mean field level with an in¬ 
termediate interaction strength 7 = 7 (see also Appendix. IbTi. 


V. DENSITY PROEILE 

In order to further test our approach, we computed the pro¬ 
file of the falling atoms at long times after they exit the ID gas, 
which can be cornpared with the experimental results reported 
in Fig. 5 of Ref. Q]. Experimental results are available |[3l for 
all three cases in Table|T] Unfortunately, a similar comparison 
was not carried out in earlier theoretical modelling lUl fTMuOr . 

In our simulations, we approximate the ID gas to have uni¬ 
form density and fixed length, with parameters Np and L as in 
the experiments. Once again, we find that the resulting den¬ 
sity profile of the falling packet has a significant sensitivity on 
the value of 7 , which allows us to readily identify which one 
gives the best fit. 

The outcome is shown in Fig.|2l The experimental and sim¬ 
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FIG. 1. (Colour online) — Position of the center of mass as a 
function of time from the numerical simulations in comparison with 
the experimental data from Fig. 3 in Ref. [U (blue open squares). The 
free-fall analytic solution a;(f) = gf ^/2 is also shown for comparison 
(red solid dots). Top panel: numerical results at fixed density n ~ 
1.2 for 7 = 1.9 (cyan open circles), 7 = 2.6 (cyan solid 

dots), and 7 = 4.0 (cyan open triangles). Bottom panel: numerical 
results at fixed 7 = 2.6, for n ~ 0.84 pm”^ (cyan open circles), 
n ~ 1.2 pm“^(cyan solid dots), and n ~ 1.56 gm~^ (cyan open 
triangles). 


ulation curves were normalised so that the area under the pro¬ 
files equals 1 (after subtraction of a background 1^ ). The 
main peak in the figure is due to the fraction of part icles that 
fall freely through the ID gas without scattering ll3ll] . 

Both the overall shape of the curves and the ratio be¬ 
tween scattered and free-falling contributions are in reason¬ 
able agreement between numerics and experiments for 7 = 
1.1, 1.9, and 2.1, respectively. These results suggest that the 
relevant values of 7 in the experiments are those from the cen¬ 
tral region of the condensate. 

We note that for these values of 7 we find very good agree¬ 
ment between the exact LL solution and the Bose gas (BG) 
approximation 1^ . In the BG limit, we studied also ID gases 
with static position-dependent density llz^ . We found that the 
resulting effects are minor and do not alter the best fit values 
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FIG. 2. Behaviour starting from gaussian distributed initial condi¬ 
tions over Ai'hist = 100, 000 histories, for 7 = 1.1, 1.9, 2.1 (black, 
red and blue, respectively) with a finite trap of uniform density (see 
text for size and density details). The results are expressed as his¬ 
tograms of the position of the particles after they have been falling 
for 18.6 ms The corresponding experimental results in Fig. 5 
of Ref. [3 are shown as thin dotted lines. 


of 7 . 

We notice that a small dip between the scattered peak and 
the free-falling peak appears in the experimental data (most 
noticeably at larger values of 7 ) whereas it is nearly absent in 
the numerical simulations. We conjecture that this dip might 
be due to the fact that in the actual experiment two or more 
impurities might fall though the trap at the same time. An 
effective attractive interaction between impurities could bind 
together nearby impurities and enhance the main peak at the 
expense of weight on either sides of the main peak. This effect 
is beyond our approximation and must be relegated to future 
research. 


VI. CONCLUSIONS 

Using linear response theory and Fermi’s golden rule with 
exact transition rates to model the scattering between the 
driven impurity and the underlying ID Bose gas, we have been 
able to obtain a quantitative description of the experimental 
results in Ref. [I]: center of mass, profile of driven packet with 
time-of-flight measurements and tomography. It appears that 
our crude approach is entirely sufficient to describe the vast 
majority of the observed physics. If desired, the approach 
taken here could be systematically improved by considering 
corrections of higher order in the coupling between the Bose 
gas and the impurity. The detection of finer quantum mechani¬ 
cal effects beyond our description may however require higher 
experimental accuracy. 

In the range of 7 values considered here, the need for an ex¬ 
act LL solution was limited and the results would have been 
in large part the same had we used the BG approximation 
instead. However, sizeable differences between LL and BG 


arise already for 7 > 3, which is within experimental reach. 

The case considered here can be viewed as an extremely 
simple example of handling a non-equilibrium situation in the 
presence of strong correlations. While the impurity is driven, 
the physics of the Bose gas can still be understood as remain¬ 
ing at equilibrium. Going further, one could try to understand 
how putting the gas itself out of equilibrium affects the im¬ 
purity dynamics. Moreover, besides cold atom settings, one 
could also consider driven quantum magnets, for which the 
necessary exact correlators are also available. We will return 
to these issues in future work. 
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Appendix A: Dependence on ID gas density profile 


Por the values of 7 of relevance to the experiments in Ref.[l], 
we find that our simulations give similar results whether we 
use the DSP from LL or in the BG approximation. We can 
therefore use the latter to test how the results are affected by 
(static) changes in the ID gas density profile. 

The DSP of a ID bose gas can be determined directly from 
its spectrum ll^ 
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After a few lines of algebra, following the steps outlined in 
Sec.ni one obtains that the only allowed outgoing wave vector 
is k' = jkp/{TT'^k), with probability 


r, 




> VtA (A3) 

0 otherwise. 


Notice that Eq. (IA3I) can be interpreted as a probability only 
if it is < 1, which in turn is satisfied if we choose 
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where we used explicitly the condition |fc|//cF > 
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Using Eq. (IA3b one can straightforwardly adapt the simu¬ 
lations to a position-dependent (static) density profile of the 
underlying ID gas. For concreteness we fix the average den¬ 
sity at the experimentally relevant value of 1.3278 particles 
per micrometre (corresponding, in the case of uniform den¬ 
sity, to an average 7 = 1.9). We then contrast the following 
cases: (i) a uniform condensate of finite length L = 24.1 /rm; 
(ii) a uniform condensate of the same length with a square de¬ 
pletion to half its density near its center (defined as — 2.0 < 
X < 2.0 /rm); (iii) a parabolic condensate of the same length 
and average density; and finally (iv) a parabolic condensate 
with a central depletion obtained by subtracting the Gaussian- 
smoothed Wigner function that we used to describe the initial 
distribution of the falling packet, Eq. Q, after setting p = 0. 
The four different options are illustrated in Fig. [3 The de- 



FIG. 3. Different (static) density profiles for the ID condensate 
used to assess how the shape affects the results of our simulations: 
(i) square (red), (ii) square with a square depletion at the center of the 
condensate (magenta), (iii) parabolic (blue), and (iv) parabolic with 
a Gaussian-smoothed Wigner function depletion at the center of the 
condensate (green). 


Appendix B: Tonks-Girardeau limit 


In our work we have found quantitative agreement with the 
experimental results in Ref. Q] for small values of 7 where the 
BG approximation is reasonably accurate. This is in contrast 
with the modelling presented in that very same reference lH, 
as well as the work done in Ref. [ij and Ref. |2^ which make 
use of the Tonks-Girardeau (TG) limit. 

In this section we investigate the motion of the center of 
mass of the falling packet in the TG limit using our method. 
A reasonable agreement with the experimental results can be 
obtained only in the small 7 limit, which is in contradiction 
with the TG approximation. According to our simulations, al¬ 
ready at intermediate values of the coupling strength (namely, 
7 > 7) the falling atoms reach terminal velocity well within 
the time of the experiment, in contrast with the observed be¬ 
haviour. 

The dynamic structure factor of a ID gas in the TG 
limit iH (7 —>■ 00 ) can be written as 


S(k,uj) = [©(w - w_)0(a;+ - w) 


— 0(w — w+)0(w_ — w)], (Bl) 


where 


UJ±{k) = 


hkF 

2m 
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Introducing the dimensionless wave vector notation k = 
k/kp, after the usual substitution k = k' — k and uj = 
h^k' / 2 m — h^k^/2m, a few lines of algebra show that the 
scattering probability density per unit of dimensionless wave 
vector, from kXok', is given by 


Vk^k' 





1 

\k'-k\ 


if Ifcl > 1 and Ifc'l < 1 
II II 


otherwise. 


pleted cases are intended to mimic the effect of the decou¬ 
pling laser that excites some of the atoms in the condensate to 
a non-trapped state, thus creating the initial packet of falling 
atoms (cf. Fig. 2 in Ref. [I]). Note that we continue to neglect 
any feedback between the falling atoms and the condensate, 
nor we allow the latter to relax to a shape different form the 
initial one. 

Fig.|4]shows the density profiles of the falling atoms at dif¬ 
ferent times, using the same initial conditions discussed in the 
main text. The differences are minor and comparable to the 
experimental error bars in Ref. Ql The case of a parabolic 
profile ought to be considered with care, since a continuously 
vanishing density at its edges implies large values of 7 , and 
the BG approximation is no longer justified. 

We notice that Ref. which considers a similar system 
in presence of a ID optical lattice, also reported qualitatively 
similar results whether the ID gas was prepared in equilib¬ 
rium with or without the impurities (see the third paragraph in 
Sec.IIDofRef.dl). 


The expression above, which is correct to leading order in 
7 , presents the intrinsic problem that the total scattering prob¬ 
ability at a given time. 


+ 00 


Vk^k' dk' 


7^ . 


In 
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diverges in the limit k ^ kp. For the stochastic approach 
to be valid, a necessary condition is that 5t be small enough 
so that the integrated probability at any given time remains 
smaller than 1 , which thus requires 5t to be vanishingly small 
for k arbitrarily close to kp. 

The singularity is directly related to the limit 7 — 00 . 
However, it cannot be easily resolved by including the sub¬ 
leading correction in 1 /j because the expansion of S becomes 
negative in some range of k and uj . 

A compromise to obtain a non-negative, non-divergent 
probability is to use the expansion of S to leading order, as 
in Eq. (IBlb . but to replace the Heaviside Theta functions with 
those from the Random Phase Approximation (RPA). Namely, 
we use the TG form of the DSF, but with support in k and oj 
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FIG. 4. Density profiles at different times starting from the initial conditions discussed in the main text, averaged over A^hist = 100, 000 
histories. From left to right: t = 1.5 ms, t = 10 ms, and t = 18.7 ms. The colours correspond to the four cases discussed in the text; (i) 
a finite uniform condensate (black); (ii) a finite uniform condensate with a square depletion at the center (red); (iii) a parabolic condensate 
(blue); and (iv) a parabolic condensate with a Gaussian-smoothed Wigner function depletion at the center (green). 


from RPA. This in turn means that the probability Vk^k' re¬ 
tains the same form as in Eq. (IB3b . but it is set to zero identi¬ 
cally outside the range: 


< or 
[k > 1 


7-4-fc(7-2) < fc' < 1 _ 2(l+fc) 

(B5) 

7-4+2fc ^ 7/ ^ 1 _ 2(l-rfc) 

^_2 K ^ L ^_2 , 


and similarly for negative values of k. In the limit of 7 —oo 
this tends to Eq. (1B31) . as one would expect. We tested our 
choice of probability regularisation in the TG limit by compar¬ 
ing its results with RPA and LL simulations for large values 
of 7 and we found good quantitative agreement (not shown). 

Using the new boundaries in Eq. (1B5I) . the probability that 
a particle with wave vector k scatters with the condensate in 
a time interval St (to any allowed k') remains finite for all 
allowed values of k. Namely, 


r-l-2(l-lfc)/(7-2) 


Vk^k' dk' 


' [7-4-fc(7-2)]/2 


TT^ h 



if 1 — 4/7 < fc < 1, and 

fc -F 1 — 4/7 
fc — 1 -F 4/7 


l_2(l+fe)/(7-2) 

/-(7-4+2fc)/(7-2) 


Vk^k' dk' = 

TT'^ a 


if fc > 1 . Notice that the maximum over fc > 1 of the logarith¬ 
mic contribution in the second case is in fact the same as the 
first case: ln [(7 — 2)/2]. Our stochastic approach is therefore 
valid, provided that we choose 


St < 


7^ EF 


In 
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Eor the typical system parametres considered in this work, the 
upper bound for St scales as ( 7 ^ 107 )“^ milliseconds. This is 
satisfied for instance by choosing St < 0.01 iJ,s up to 7 = 100. 

We can then implement our stochastic approach using the 
inverse transform sampling analytically in the TG limit. Eig.|5] 
shows the resulting behaviour of the centre of mass (CM) mo¬ 
tion from TG simulations for a uniform ID gas of density 


n = 1.2 /im ^ and 7 = 4.1, 5, 6, 8, 10, 25, 100, to be con¬ 
trasted with the results presented in the main text, Eig.[T] 

We notice that the CM motion becomes asymptotically lin¬ 
ear in time within the simulation time window for 7 > 7, 
suggesting that the falling atoms reach terminal velocity. The 
value of the terminal velocity is non-monotonic in 7 : it ini¬ 
tially decreases (in agreement with Ref. [fsl) with increasing 
7 , and later increases and tends asymptotically to vp in the 
7 —>■ 00 limit, as expected. 

Reasonable agreement with the experimental results can 
only be obtained in the weak coupling limit (7 ~ 5), which is 
in contradiction with the TG limit (and even with the RPA ap¬ 
proximation, which has a hard limit of applicability of 7 > 4, 
and is known to begin to fit reasonably well the LL DSE only 

for7 > 10 iH). 



FIG. 5. (Colour online) — Position of the center of mass as a func¬ 
tion of time from our simulations in the TG limit, considering a uni¬ 
form ID gas with n = 1.2 /rm“^ and 7 = 4.1, 5, 6 , 8 , 10, 25, 100 
(red, blue, green, magenta, cyan, yellow, and black, respectively). 
The black dashed line represents the free-fall curve. The black solid 
line corresponds to the expected behaviour in the 7 —>■ 00 limit (i.e., 
terminal velocity vf)- Blue open squares (joined^ a dotted line) 
represent the experimental data from Fig. 3 in Ref.lJ. 
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